home *** CD-ROM | disk | FTP | other *** search
/ Aminet 35 / Aminet 35 (2000)(Schatztruhe)[!][Feb 2000].iso / Aminet / gfx / misc / gnuplot-src.lha / gnuplot-3.7.1src / gnuplot-3.7.1.lha / gnuplot-3.7.1 / interpol.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-08-19  |  30.4 KB  |  1,044 lines
  1. #ifndef lint    } elSid = "$Id: interpol.c,v 1.6.2.1 1999/08/19 14:39:18 lhecking Exp $";
  2. #endif
  3.  
  4. /* GNUPLOT - interpol.c */
  5.  
  6. /*[
  7.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  8.  *
  9.  * Permission to use, copy, and distribute this software and its
  10.  * documentation for any purpose with or without fee is hereby granted,
  11.  * provided that the above copyright notice appear in all copies and
  12.  * that both that copyright notice and this permission notice appear
  13.  * in supporting documentation.
  14.  *
  15.  * Permission to modify the software is granted, but not the right to
  16.  * distribute the complete modified source code.  Modifications are to
  17.  * be distributed as patches to the released version.  Permission to
  18.  * distribute binaries produced by compiling modified sources is granted,
  19.  * provided you
  20.  *   1. distribute the corresponding source modifications from the
  21.  *    released version in the form of a patch file along with the binaries,
  22.  *   2. add special version identification to distinguish your version
  23.  *    in addition to the base release version number,
  24.  *   3. provide your name and address as the primary contact for the
  25.  *    support of your modified version, and
  26.  *   4. retain our contact information in regard to use of the base
  27.  *    software.
  28.  * Permission to distribute the released version of the source code along
  29.  * with corresponding source modifications in the form of a patch file is
  30.  * granted with same provisions 2 through 4 for binary distributions.
  31.  *
  32.  * This software is provided "as is" without express or implied warranty
  33.  * to the extent permitted by applicable law.
  34. ]*/
  35.  
  36.  
  37. /*
  38.  * C-Source file identification Header
  39.  *
  40.  * This file belongs to a project which is:
  41.  *
  42.  * done 1993 by MGR-Software, Asgard  (Lars Hanke)
  43.  * written by Lars Hanke
  44.  *
  45.  * Contact me via:
  46.  *
  47.  *  InterNet: mgr@asgard.bo.open.de
  48.  *      FIDO: Lars Hanke @ 2:243/4802.22   (as long as they keep addresses)
  49.  *
  50.  **************************************************************************
  51.  *
  52.  *   Project: gnuplot
  53.  *    Module:
  54.  *      File: interpol.c
  55.  *
  56.  *   Revisor: Lars Hanke
  57.  *   Revised: 26/09/93
  58.  *  Revision: 1.0
  59.  *
  60.  **************************************************************************
  61.  *
  62.  * LEGAL
  63.  *  This module is part of gnuplot and distributed under whatever terms
  64.  *  gnuplot is or will be published, unless exclusive rights are claimed.
  65.  *
  66.  * DESCRIPTION
  67.  *  Supplies 2-D data interpolation and approximation routines
  68.  *
  69.  * IMPORTS
  70.  *  plot.h
  71.  *    - cp_extend()
  72.  *    - structs: curve_points, coordval, coordinate
  73.  *
  74.  *  setshow.h
  75.  *    - samples, is_log_x, base_log_x, xmin, xmax, autoscale_lx
  76.  *    - plottypes
  77.  *
  78.  *  proto.h
  79.  *    - solve_tri_diag()
  80.  *    - typedef tri_diag
  81.  *
  82.  * EXPORTS
  83.  *  gen_interp()
  84.  *  sort_points()
  85.  *  cp_implode()
  86.  *
  87.  * BUGS and TODO
  88.  *  I would really have liked to use Gershon Elbers contouring code for
  89.  *  all the stuff done here, but I failed. So I used my own code.
  90.  *  If somebody is able to consolidate Gershon's code for this purpose
  91.  *  a lot of gnuplot users would be very happy - due to memory problems.
  92.  *
  93.  **************************************************************************
  94.  *
  95.  * HISTORY
  96.  * Changes:
  97.  *  Nov 24, 1995  Markus Schuh (M.Schuh@meteo.uni-koeln.de):
  98.  *      changed the algorithm for csplines
  99.  *      added algorithm for approximation csplines
  100.  *      copied point storage and range fix from plot2d.c
  101.  *
  102.  *  Dec 12, 1995 David Denholm
  103.  *      oops - at the time this is called, stored co-ords are
  104.  *      internal (ie maybe log of data) but min/max are in
  105.  *      user co-ordinates. 
  106.  *      Work with min and max of internal co-ords, and
  107.  *      check at the end whether external min and max need to
  108.  *      be increased. (since samples is typically 100 ; we
  109.  *      dont want to take more logs than necessary)
  110.  *      Also, need to take into account which axes are active
  111.  *
  112.  *  Jun 30, 1996 Jens Emmerich
  113.  *      implemented handling of UNDEFINED points
  114.  */
  115.  
  116. #include "plot.h"
  117. #include "setshow.h"
  118.  
  119.  
  120. /* in order to support multiple axes, and to simplify ranging in
  121.  * parametric plots, we use arrays to store some things. For 2d plots,
  122.  * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
  123.  * names in plot.h
  124.  */
  125.  
  126. extern double min_array[AXIS_ARRAY_SIZE];
  127. extern double max_array[AXIS_ARRAY_SIZE];
  128. extern int auto_array[AXIS_ARRAY_SIZE];
  129. extern TBOOLEAN log_array[AXIS_ARRAY_SIZE];
  130. extern double base_array[AXIS_ARRAY_SIZE];
  131. extern double log_base_array[AXIS_ARRAY_SIZE];
  132.  
  133.  
  134. #define Inc_c_token if (++c_token >= num_tokens)    \
  135. int_error ("Syntax error", c_token);
  136.  
  137.  
  138. /*
  139.  * IMHO, code is getting too cluttered with repeated chunks of
  140.  * code. Some macros to simplify, I hope.
  141.  *
  142.  * do { } while(0) is comp.lang.c recommendation for complex macros
  143.  * also means that break can be specified as an action, and it will
  144.  * 
  145.  */
  146.  
  147.  
  148. /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
  149.  * Do OUT_ACTION or UNDEF_ACTION as appropriate
  150.  * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
  151.  * VALUE must not be same as STORE
  152.  */
  153.  
  154. #define STORE_AND_FIXUP_RANGE(STORE, VALUE, TYPE, MIN, MAX, AUTO, OUT_ACTION, UNDEF_ACTION)\
  155. do { STORE=VALUE; \
  156.     if (TYPE != INRANGE) break;  /* dont set y range if x is outrange, for example */ \
  157.     if ( VALUE<MIN ) { \
  158.        if (AUTO & 1) MIN=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; break; }  \
  159.     } \
  160.     if ( VALUE>MAX ) {\
  161.        if (AUTO & 2) MAX=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; }   \
  162.     } \
  163. } while(0)
  164.  
  165. #define UPDATE_RANGE(TEST,OLD,NEW,AXIS) \
  166. do { if (TEST) { \
  167.      if (log_array[AXIS]) OLD = pow(base_array[AXIS], NEW); else OLD = NEW; \
  168.      } \
  169. } while(0)
  170.  
  171. /* use this instead empty macro arguments to work around NeXT cpp bug */
  172. /* if this fails on any system, we might use ((void)0) */
  173.  
  174. #define NOOP            /* */
  175.  
  176. #define spline_coeff_size 4
  177. typedef double spline_coeff[spline_coeff_size];
  178. typedef double five_diag[5];
  179.  
  180. static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
  181. static int num_curves __PROTO((struct curve_points * plot));
  182. static double *cp_binomial __PROTO((int points));
  183. GP_INLINE static double s_pow __PROTO((double base, unsigned int exponent));
  184. static void eval_bezier __PROTO((struct curve_points * cp, int first_point, int num_points, double sr, coordval * px, coordval * py, double *c));
  185. static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
  186. static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
  187. static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
  188. static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
  189. static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
  190. static int compare_points __PROTO((struct coordinate * p1, struct coordinate * p2));
  191.  
  192.  
  193. /*
  194.  * position curve_start to index the next non-UNDEFINDED point,
  195.  * start search at initial curve_start,
  196.  * return number of non-UNDEFINDED points from there on,
  197.  * if no more valid points are found, curve_start is set
  198.  * to plot->p_count and 0 is returned
  199.  */
  200.  
  201. static int next_curve(plot, curve_start)
  202. struct curve_points *plot;
  203. int *curve_start;
  204. {
  205.     int curve_length;
  206.  
  207.     /* Skip undefined points */
  208.     while (*curve_start < plot->p_count
  209.        && plot->points[*curve_start].type == UNDEFINED) {
  210.     (*curve_start)++;
  211.     };
  212.     curve_length = 0;
  213.     /* curve_length is first used as an offset, then the correkt # points */
  214.     while ((*curve_start) + curve_length < plot->p_count
  215.       && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
  216.     curve_length++;
  217.     };
  218.     return (curve_length);
  219. }
  220.  
  221.  
  222. /* 
  223.  * determine the number of curves in plot->points, separated by
  224.  * UNDEFINED points
  225.  */
  226.  
  227. static int num_curves(plot)
  228. struct curve_points *plot;
  229. {
  230.     int curves;
  231.     int first_point;
  232.     int num_points;
  233.  
  234.     first_point = 0;
  235.     curves = 0;
  236.     while ((num_points = next_curve(plot, &first_point)) > 0) {
  237.     curves++;
  238.     first_point += num_points;
  239.     }
  240.     return (curves);
  241. }
  242.  
  243.  
  244.  
  245. /*
  246.  * build up a cntr_struct list from curve_points
  247.  * this funtion is only used for the alternate entry point to
  248.  * Gershon's code and thus commented out
  249.  ***deleted***
  250.  */
  251.  
  252.  
  253. /* HBB 990205: rewrote the 'bezier' interpolation routine,
  254.  * to prevent numerical overflow and other undesirable things happening
  255.  * for large data files (num_data about 1000 or so), where binomial
  256.  * coefficients would explode, and powers of 'sr' (0 < sr < 1) become
  257.  * extremely small. Method used: compute logarithms of these
  258.  * extremely large and small numbers, and only go back to the
  259.  * real numbers once they've cancelled out each other, leaving
  260.  * a reasonable-sized one. */
  261.  
  262. /*
  263.  * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
  264.  *   and stores them into an array which is hooked to sdat.
  265.  * (MGR 1992)
  266.  */
  267. static double *cp_binomial(points)
  268. int points;
  269. {
  270.     register double *coeff;
  271.     register int n, k;
  272.     int e;
  273.  
  274.     e = points;            /* well we're going from k=0 to k=p_count-1 */
  275.     coeff = (double *) gp_alloc(e * sizeof(double), "bezier coefficients");
  276.  
  277.     n = points - 1;
  278.     e = n / 2;
  279.     /* HBB 990205: calculate these in 'logarithmic space',
  280.      * as they become _very_ large, with growing n (4^n) */
  281.     coeff[0] = 0.0;
  282.  
  283.     for (k = 0; k < e; k++) {
  284.     coeff[k + 1] = coeff[k] + log(((double) (n-k)) / ((double) (k+1)));
  285.     }
  286.  
  287.     for (k = n; k >= e; k--)
  288.     coeff[k] = coeff[n - k];
  289.  
  290.     return (coeff);
  291. }
  292.  
  293.  
  294. /* This is a subfunction of do_bezier() for BEZIER style computations.
  295.  * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
  296.  * the double values holding the next x and y coordinates.
  297.  * (MGR 1992)
  298.  */
  299.  
  300. /*
  301.  * well, this routine runs faster with the 68040 striptease FPU
  302.  * and it handles zeroes correctly - there had been some trouble with TC
  303.  */
  304.  
  305. GP_INLINE static double s_pow(base, exponent)
  306. double base;
  307. unsigned int exponent;
  308. {
  309.     double y;
  310.  
  311.     if (exponent == 0)
  312.     return (1.0);
  313.     if (base == 0.0)
  314.     return (0.0);
  315.  
  316.     /* consider i in binary = abcd
  317.      * x^i = x^(8a+4b+2c+d) = x^8a * x^4b * x^2b * x^d
  318.      *                      = a?x^2^2^2:1 * b?x^2^2:1 + ...
  319.      * so for each bit in exponent, square x, multiplying it into accumulator
  320.      *
  321.      */
  322.  
  323.     y = 1;
  324.     while (exponent) {
  325.     if (exponent & 1)
  326.         y *= base;
  327.     base *= base;
  328.     /* if exponent was signed, this could be trouble ! */
  329.     exponent >>= 1;
  330.     }
  331.  
  332.     return (y);
  333. }
  334.  
  335.  
  336. static void eval_bezier(cp, first_point, num_points, sr, px, py, c)
  337. struct curve_points *cp;
  338. int first_point;        /* where to start in plot->points (to find x-range) */
  339. int num_points;            /* to determine end in plot->points */
  340. double sr;
  341. coordval *px;
  342. coordval *py;
  343. double *c;
  344. {
  345.     unsigned int n = num_points - 1;
  346.     /* HBB 980308: added 'GPHUGE' tag for DOS */
  347.     struct coordinate GPHUGE *this_points;
  348.  
  349.     this_points = (cp->points) + first_point;
  350.  
  351.     if (sr == 0.0) {
  352.     *px = this_points[0].x;
  353.     *py = this_points[0].y;
  354.     } else if (sr == 1.0) {
  355.     *px = this_points[n].x;
  356.     *py = this_points[n].y;
  357.     } else {
  358.         /* HBB 990205: do calculation in 'logarithmic space',
  359.          * to avoid over/underflow errors, which would exactly cancel
  360.          * out each other, anyway, in an exact calculation
  361.          */
  362.     unsigned int i;
  363.     double lx = 0.0, ly = 0.0;
  364.     double log_dsr_to_the_n = n * log(1 -sr);
  365.     double log_sr_over_dsr = log(sr) - log(1 - sr);
  366.  
  367.     for (i = 0; i <= n; i++) {
  368.             double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr);
  369.  
  370.         lx += this_points[i].x * u;
  371.         ly += this_points[i].y * u;
  372.     }
  373.  
  374.     *px = lx;
  375.     *py = ly;
  376.     }
  377. }
  378.  
  379. /*
  380.  * generate a new set of coordinates representing the bezier curve and
  381.  * set it to the plot
  382.  */
  383.  
  384. static void do_bezier(cp, bc, first_point, num_points, dest)
  385. struct curve_points *cp;
  386. double *bc;
  387. int first_point;        /* where to start in plot->points */
  388. int num_points;            /* to determine end in plot->points */
  389. struct coordinate *dest;    /* where to put the interpolated data */
  390. {
  391.     int i;
  392.     coordval x, y;
  393.     int xaxis = cp->x_axis;
  394.     int yaxis = cp->y_axis;
  395.  
  396.     /* min and max in internal (eg logged) co-ordinates. We update
  397.      * these, then update the external extrema in user co-ordinates
  398.      * at the end.
  399.      */
  400.  
  401.     double ixmin, ixmax, iymin, iymax;
  402.     double sxmin, sxmax, symin, symax;    /* starting values of above */
  403.  
  404.     if (log_array[xaxis]) {
  405.     ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
  406.     ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
  407.     } else {
  408.     ixmin = sxmin = min_array[xaxis];
  409.     ixmax = sxmax = max_array[xaxis];
  410.     }
  411.  
  412.     if (log_array[yaxis]) {
  413.     iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
  414.     iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
  415.     } else {
  416.     iymin = symin = min_array[yaxis];
  417.     iymax = symax = max_array[yaxis];
  418.     }
  419.  
  420.     for (i = 0; i < samples; i++) {
  421.     eval_bezier(cp, first_point, num_points, (double) i / (double) (samples - 1), &x, &y, bc);
  422.  
  423.     /* now we have to store the points and adjust the ranges */
  424.  
  425.     dest[i].type = INRANGE;
  426.     STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
  427.     STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
  428.  
  429.     dest[i].xlow = dest[i].xhigh = dest[i].x;
  430.     dest[i].ylow = dest[i].yhigh = dest[i].y;
  431.  
  432.     dest[i].z = -1;
  433.     }
  434.  
  435.     UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
  436.     UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
  437.     UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
  438.     UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
  439. }
  440.  
  441. /*
  442.  * call contouring routines -- main entry
  443.  */
  444.  
  445. /* 
  446.  * it should be like this, but it doesn't run. If you find out why,
  447.  * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
  448.  *
  449.  * Well, all this had originally been inside contour.c, so maybe links
  450.  * to functions and of contour.c are broken.
  451.  * ***deleted***
  452.  * end of unused entry point to Gershon's code
  453.  *
  454.  */
  455.  
  456. /* 
  457.  * Solve five diagonal linear system equation. The five diagonal matrix is
  458.  * defined via matrix M, right side is r, and solution X i.e. M * X = R.
  459.  * Size of system given in n. Return TRUE if solution exist.
  460.  *  G. Engeln-Muellges/ F.Reutter: 
  461.  *  "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
  462.  *  ISBN 3-411-01677-9
  463.  *
  464.  * /  m02 m03 m04   0   0   0   0    .       .       . \   /  x0  \    / r0  \
  465.  * I  m11 m12 m13 m14   0   0   0    .       .       . I   I  x1  I   I  r1  I
  466.  * I  m20 m21 m22 m23 m24   0   0    .       .       . I * I  x2  I = I  r2  I
  467.  * I    0 m30 m31 m32 m33 m34   0    .       .       . I   I  x3  I   I  r3  I
  468.  *      .   .   .   .   .   .   .    .       .       .        .        .
  469.  * \                           m(n-3)0 m(n-2)1 m(n-1)2 /   \x(n-1)/   \r(n-1)/
  470.  * 
  471.  */
  472. static int solve_five_diag(m, r, x, n)
  473. five_diag m[];
  474. double r[], x[];
  475. int n;
  476. {
  477.     int i;
  478.     five_diag *hv;
  479.  
  480.     hv = (five_diag *) gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
  481.  
  482.     hv[0][0] = m[0][2];
  483.     if (hv[0][0] == 0) {
  484.     free(hv);
  485.     return FALSE;
  486.     }
  487.     hv[0][1] = m[0][3] / hv[0][0];
  488.     hv[0][2] = m[0][4] / hv[0][0];
  489.  
  490.     hv[1][3] = m[1][1];
  491.     hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
  492.     if (hv[1][0] == 0) {
  493.     free(hv);
  494.     return FALSE;
  495.     }
  496.     hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
  497.     hv[1][2] = m[1][4] / hv[1][0];
  498.  
  499.     for (i = 2; i <= n - 1; i++) {
  500.     hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
  501.     hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
  502.     if (hv[i][0] == 0) {
  503.         free(hv);
  504.         return FALSE;
  505.     }
  506.     hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
  507.     hv[i][2] = m[i][4] / hv[i][0];
  508.     }
  509.  
  510.     hv[0][4] = 0;
  511.     hv[1][4] = r[0] / hv[0][0];
  512.     for (i = 1; i <= n - 1; i++) {
  513.     hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
  514.     }
  515.  
  516.     x[n - 1] = hv[n][4];
  517.     x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
  518.     for (i = n - 3; i >= 0; i--)
  519.     x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
  520.  
  521.     free(hv);
  522.     return TRUE;
  523. }
  524.  
  525.  
  526. /*
  527.  * Calculation of approximation cubic splines
  528.  * Input:  x[i], y[i], weights z[i]
  529.  *         
  530.  * Returns matrix of spline coefficients
  531.  */
  532. static spline_coeff *cp_approx_spline(plot, first_point, num_points)
  533. struct curve_points *plot;
  534. int first_point;        /* where to start in plot->points */
  535. int num_points;            /* to determine end in plot->points */
  536. {
  537.     spline_coeff *sc;
  538.     five_diag *m;
  539.     int xaxis = plot->x_axis;
  540.     int yaxis = plot->y_axis;
  541.     double *r, *x, *h, *xp, *yp;
  542.     /* HBB 980308: added 'GPHUGE' tag */
  543.     struct coordinate GPHUGE *this_points;
  544.     int i;
  545.  
  546.     sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff),
  547.                    "spline matrix");
  548.  
  549.     if (num_points < 4)
  550.     int_error("Can't calculate approximation splines, need at least 4 points", NO_CARET);
  551.  
  552.     this_points = (plot->points) + first_point;
  553.  
  554.     for (i = 0; i <= num_points - 1; i++)
  555.     if (this_points[i].z <= 0)
  556.         int_error("Can't calculate approximation splines, all weights have to be > 0", NO_CARET);
  557.  
  558.     m = (five_diag *) gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
  559.  
  560.     r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
  561.     x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
  562.     h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
  563.  
  564.     xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
  565.     yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
  566.  
  567.     /* KB 981107: With logarithmic axis first convert back to linear scale */
  568.  
  569.     if (log_array[xaxis]) {
  570.     for (i = 0; i <= num_points - 1; i++)
  571.         xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
  572.     } else {
  573.     for (i = 0; i <= num_points - 1; i++)
  574.         xp[i] = this_points[i].x;
  575.     }
  576.     if (log_array[yaxis]) {
  577.     for (i = 0; i <= num_points - 1; i++)
  578.         yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
  579.     } else {
  580.     for (i = 0; i <= num_points - 1; i++)
  581.         yp[i] = this_points[i].y;
  582.     }
  583.  
  584.     for (i = 0; i <= num_points - 2; i++)
  585.     h[i] = xp[i + 1] - xp[i];
  586.  
  587.     /* set up the matrix and the vector */
  588.  
  589.     for (i = 0; i <= num_points - 3; i++) {
  590.     r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
  591.             - (yp[i + 1] - yp[i]) / h[i]);
  592.  
  593.     if (i < 2)
  594.         m[i][0] = 0;
  595.     else
  596.         m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
  597.  
  598.     if (i < 1)
  599.         m[i][1] = 0;
  600.     else
  601.         m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
  602.         - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
  603.  
  604.     m[i][2] = 2 * (h[i] + h[i + 1])
  605.         + 6 / this_points[i].z / h[i] / h[i]
  606.         + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
  607.         + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
  608.  
  609.     if (i > num_points - 4)
  610.         m[i][3] = 0;
  611.     else
  612.         m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
  613.         - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
  614.  
  615.     if (i > num_points - 5)
  616.         m[i][4] = 0;
  617.     else
  618.         m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
  619.     }
  620.  
  621.     /* solve the matrix */
  622.     if (!solve_five_diag(m, r, x, num_points - 2)) {
  623.     free(h);
  624.     free(x);
  625.     free(r);
  626.     free(m);
  627.     free(xp);
  628.     free(yp);
  629.     int_error("Can't calculate approximation splines", NO_CARET);
  630.     }
  631.     sc[0][2] = 0;
  632.     for (i = 1; i <= num_points - 2; i++)
  633.     sc[i][2] = x[i - 1];
  634.     sc[num_points - 1][2] = 0;
  635.  
  636.     sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
  637.     for (i = 1; i <= num_points - 2; i++)
  638.     sc[i][0] = yp[i] - 2 / this_points[i].z *
  639.         (sc[i - 1][2] / h[i - 1]
  640.          - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
  641.          + sc[i + 1][2] / h[i]);
  642.     sc[num_points - 1][0] = yp[num_points - 1]
  643.     - 2 / this_points[num_points - 1].z / h[num_points - 2]
  644.     * (sc[num_points - 2][2] - sc[num_points - 1][2]);
  645.  
  646.     for (i = 0; i <= num_points - 2; i++) {
  647.     sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
  648.         - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
  649.     sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
  650.     }
  651.  
  652.     free(h);
  653.     free(x);
  654.     free(r);
  655.     free(m);
  656.     free(xp);
  657.     free(yp);
  658.  
  659.     return (sc);
  660. }
  661.  
  662.  
  663. /*
  664.  * Calculation of cubic splines
  665.  *
  666.  * This can be treated as a special case of approximation cubic splines, with
  667.  * all weights -> infinity.
  668.  *         
  669.  * Returns matrix of spline coefficients
  670.  */
  671. static spline_coeff *cp_tridiag(plot, first_point, num_points)
  672. struct curve_points *plot;
  673. int first_point, num_points;
  674. {
  675.     spline_coeff *sc;
  676.     tri_diag *m;
  677.     int xaxis = plot->x_axis;
  678.     int yaxis = plot->y_axis;
  679.     double *r, *x, *h, *xp, *yp;
  680.     /* HBB 980308: added 'GPHUGE' tag */
  681.     struct coordinate GPHUGE *this_points;
  682.     int i;
  683.  
  684.     if (num_points < 3)
  685.     int_error("Can't calculate splines, need at least 3 points", NO_CARET);
  686.  
  687.     this_points = (plot->points) + first_point;
  688.  
  689.     sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
  690.     m = (tri_diag *) gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
  691.  
  692.     r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
  693.     x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
  694.     h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
  695.  
  696.     xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
  697.     yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
  698.  
  699.     /* KB 981107: With logarithmic axis first convert back to linear scale */
  700.  
  701.     if (log_array[xaxis]) {
  702.     for (i = 0; i <= num_points - 1; i++)
  703.         xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
  704.     } else {
  705.     for (i = 0; i <= num_points - 1; i++)
  706.         xp[i] = this_points[i].x;
  707.     }
  708.     if (log_array[yaxis]) {
  709.     for (i = 0; i <= num_points - 1; i++)
  710.         yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
  711.     } else {
  712.     for (i = 0; i <= num_points - 1; i++)
  713.         yp[i] = this_points[i].y;
  714.     }
  715.  
  716.     for (i = 0; i <= num_points - 2; i++)
  717.     h[i] = xp[i + 1] - xp[i];
  718.  
  719.     /* set up the matrix and the vector */
  720.  
  721.     for (i = 0; i <= num_points - 3; i++) {
  722.     r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
  723.             - (yp[i + 1] - yp[i]) / h[i]);
  724.  
  725.     if (i < 1)
  726.         m[i][0] = 0;
  727.     else
  728.         m[i][0] = h[i];
  729.  
  730.     m[i][1] = 2 * (h[i] + h[i + 1]);
  731.  
  732.     if (i > num_points - 4)
  733.         m[i][2] = 0;
  734.     else
  735.         m[i][2] = h[i + 1];
  736.     }
  737.  
  738.     /* solve the matrix */
  739.     if (!solve_tri_diag(m, r, x, num_points - 2)) {
  740.     free(h);
  741.     free(x);
  742.     free(r);
  743.     free(m);
  744.     free(xp);
  745.     free(yp);
  746.     int_error("Can't calculate cubic splines", NO_CARET);
  747.     }
  748.     sc[0][2] = 0;
  749.     for (i = 1; i <= num_points - 2; i++)
  750.     sc[i][2] = x[i - 1];
  751.     sc[num_points - 1][2] = 0;
  752.  
  753.     for (i = 0; i <= num_points - 1; i++)
  754.     sc[i][0] = yp[i];
  755.  
  756.     for (i = 0; i <= num_points - 2; i++) {
  757.     sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
  758.         - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
  759.     sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
  760.     }
  761.  
  762.     free(h);
  763.     free(x);
  764.     free(r);
  765.     free(m);
  766.     free(xp);
  767.     free(yp);
  768.  
  769.     return (sc);
  770. }
  771.  
  772. static void do_cubic(plot, sc, first_point, num_points, dest)
  773. struct curve_points *plot;    /* still containes old plot->points */
  774. spline_coeff *sc;        /* generated by cp_tridiag */
  775. int first_point;        /* where to start in plot->points */
  776. int num_points;            /* to determine end in plot->points */
  777. struct coordinate *dest;    /* where to put the interpolated data */
  778. {
  779.     double xdiff, temp, x, y;
  780.     int i, l;
  781.     int xaxis = plot->x_axis;
  782.     int yaxis = plot->y_axis;
  783.     /* HBB 980308: added 'GPHUGE' tag */
  784.     struct coordinate GPHUGE *this_points;
  785.  
  786.     /* min and max in internal (eg logged) co-ordinates. We update
  787.      * these, then update the external extrema in user co-ordinates
  788.      * at the end.
  789.      */
  790.  
  791.     double ixmin, ixmax, iymin, iymax;
  792.     double sxmin, sxmax, symin, symax;    /* starting values of above */
  793.  
  794.     if (log_array[xaxis]) {
  795.     ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
  796.     ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
  797.     } else {
  798.     ixmin = sxmin = min_array[xaxis];
  799.     ixmax = sxmax = max_array[xaxis];
  800.     }
  801.  
  802.     if (log_array[yaxis]) {
  803.     iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
  804.     iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
  805.     } else {
  806.     iymin = symin = min_array[yaxis];
  807.     iymax = symax = max_array[yaxis];
  808.     }
  809.  
  810.  
  811.     this_points = (plot->points) + first_point;
  812.  
  813.     l = 0;
  814.  
  815.     xdiff = (this_points[num_points - 1].x - this_points[0].x) / (samples - 1);
  816.  
  817.     for (i = 0; i < samples; i++) {
  818.     x = this_points[0].x + i * xdiff;
  819.  
  820.     while ((x >= this_points[l + 1].x) && (l < num_points - 2))
  821.         l++;
  822.  
  823.     /* KB 981107: With logarithmic x axis the values were converted back to linear  */
  824.     /* scale before calculating the coefficients. Use exponential for log x values. */
  825.  
  826.     if (log_array[xaxis]) {
  827.         temp = exp(x * log_base_array[xaxis]) - exp(this_points[l].x * log_base_array[xaxis]);
  828.         y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
  829.     } else {
  830.         temp = x - this_points[l].x;
  831.         y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
  832.     }
  833.     /* With logarithmic y axis, we need to convert from linear to log scale now. */
  834.     if (log_array[yaxis]) {
  835.         if (y > 0.)
  836.         y = log(y) / log_base_array[yaxis];
  837.         else
  838.         y = symin - (symax - symin);
  839.     }
  840.     dest[i].type = INRANGE;
  841.     STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
  842.     STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
  843.  
  844.     dest[i].xlow = dest[i].xhigh = dest[i].x;
  845.     dest[i].ylow = dest[i].yhigh = dest[i].y;
  846.  
  847.     dest[i].z = -1;
  848.  
  849.     }
  850.  
  851.     UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
  852.     UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
  853.     UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
  854.     UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
  855.  
  856. }
  857.  
  858. /*
  859.  * This is the main entry point used. As stated in the header, it is fine,
  860.  * but I'm not too happy with it.
  861.  */
  862.  
  863. void gen_interp(plot)
  864. struct curve_points *plot;
  865. {
  866.  
  867.     spline_coeff *sc;
  868.     double *bc;
  869.     struct coordinate *new_points;
  870.     int i, curves;
  871.     int first_point, num_points;
  872.  
  873.     curves = num_curves(plot);
  874.     new_points = (struct coordinate *) gp_alloc((samples + 1) * curves * sizeof(struct coordinate), "interpolation table");
  875.  
  876.     first_point = 0;
  877.     for (i = 0; i < curves; i++) {
  878.     num_points = next_curve(plot, &first_point);
  879.     switch (plot->plot_smooth) {
  880.     case CSPLINES:
  881.         sc = cp_tridiag(plot, first_point, num_points);
  882.         do_cubic(plot, sc, first_point, num_points,
  883.              new_points + i * (samples + 1));
  884.         free(sc);
  885.         break;
  886.     case ACSPLINES:
  887.         sc = cp_approx_spline(plot, first_point, num_points);
  888.         do_cubic(plot, sc, first_point, num_points,
  889.              new_points + i * (samples + 1));
  890.         free(sc);
  891.         break;
  892.  
  893.     case BEZIER:
  894.     case SBEZIER:
  895.         bc = cp_binomial(num_points);
  896.         do_bezier(plot, bc, first_point, num_points,
  897.               new_points + i * (samples + 1));
  898.         free((char *) bc);
  899.         break;
  900.     default:        /* keep gcc -Wall quiet */
  901.         ;
  902.     }
  903.     new_points[(i + 1) * (samples + 1) - 1].type = UNDEFINED;
  904.     first_point += num_points;
  905.     }
  906.  
  907.     free(plot->points);
  908.     plot->points = new_points;
  909.     plot->p_max = curves * (samples + 1);
  910.     plot->p_count = plot->p_max - 1;
  911.  
  912.     return;
  913. }
  914.  
  915. /* 
  916.  * sort_points
  917.  *
  918.  * sort data succession for further evaluation by plot_splines, etc.
  919.  * This routine is mainly introduced for compilers *NOT* supporting the
  920.  * UNIX qsort() routine. You can then easily replace it by more convenient
  921.  * stuff for your compiler.
  922.  * (MGR 1992)
  923.  */
  924.  
  925. static int compare_points(p1, p2)
  926. struct coordinate *p1;
  927. struct coordinate *p2;
  928. {
  929.     if (p1->x > p2->x)
  930.     return (1);
  931.     if (p1->x < p2->x)
  932.     return (-1);
  933.     return (0);
  934. }
  935.  
  936. void sort_points(plot)
  937. struct curve_points *plot;
  938. {
  939.     int first_point, num_points;
  940.  
  941.     first_point = 0;
  942.     while ((num_points = next_curve(plot, &first_point)) > 0) {
  943.     /* Sort this set of points, does qsort handle 1 point correctly? */
  944.     qsort((char *) (plot->points + first_point), num_points,
  945.           sizeof(struct coordinate), (sortfunc) compare_points);
  946.     first_point += num_points;
  947.     }
  948.     return;
  949. }
  950.  
  951.  
  952. /*
  953.  * cp_implode() if averaging is selected this function computes the new
  954.  *              entries and shortens the whole thing to the necessary
  955.  *              size
  956.  * MGR Addendum
  957.  */
  958.  
  959. void cp_implode(cp)
  960. struct curve_points *cp;
  961. {
  962.     int first_point, num_points;
  963.     int i, j, k;
  964.     double x, y, sux, slx, suy, sly;
  965.     enum coord_type dot;
  966.  
  967.  
  968.     j = 0;
  969.     first_point = 0;
  970.     while ((num_points = next_curve(cp, &first_point)) > 0) {
  971.     k = 0;
  972.     for (i = first_point; i < first_point + num_points; i++) {
  973.         if (!k) {
  974.         x = cp->points[i].x;
  975.         y = cp->points[i].y;
  976.         sux = cp->points[i].xhigh;
  977.         slx = cp->points[i].xlow;
  978.         suy = cp->points[i].yhigh;
  979.         sly = cp->points[i].ylow;
  980.         dot = INRANGE;
  981.         if (cp->points[i].type != INRANGE)
  982.             dot = UNDEFINED;    /* This means somthing other than usual *//* just signal to check if INRANGE */
  983.         k = 1;
  984.         } else if (cp->points[i].x == x) {
  985.         y += cp->points[i].y;
  986.         sux += cp->points[i].xhigh;
  987.         slx += cp->points[i].xlow;
  988.         suy += cp->points[i].yhigh;
  989.         sly += cp->points[i].ylow;
  990.         if (cp->points[i].type != INRANGE)
  991.             dot = UNDEFINED;
  992.         k++;
  993.         } else {
  994.         cp->points[j].x = x;
  995.         cp->points[j].y = y / (double) k;
  996.         cp->points[j].xhigh = sux / (double) k;
  997.         cp->points[j].xlow = slx / (double) k;
  998.         cp->points[j].yhigh = suy / (double) k;
  999.         cp->points[j].ylow = sly / (double) k;
  1000.         cp->points[j].type = INRANGE;
  1001.         if (dot != INRANGE) {
  1002.             if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
  1003.             cp->points[j].type = OUTRANGE;
  1004.             else if (!autoscale_ly) {
  1005.             if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
  1006.                 cp->points[j].type = OUTRANGE;
  1007.             } else
  1008.             cp->points[j].type = INRANGE;
  1009.         }
  1010.         j++;        /* next valid entry */
  1011.         k = 0;        /* to read */
  1012.         i--;        /* from this (-> last after for(;;)) entry */
  1013.         }
  1014.     }
  1015.     if (k) {
  1016.         cp->points[j].x = x;
  1017.         cp->points[j].y = y / (double) k;
  1018.         cp->points[j].xhigh = sux / (double) k;
  1019.         cp->points[j].xlow = slx / (double) k;
  1020.         cp->points[j].yhigh = suy / (double) k;
  1021.         cp->points[j].ylow = sly / (double) k;
  1022.         cp->points[j].type = INRANGE;
  1023.         if (dot != INRANGE) {
  1024.         if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
  1025.             cp->points[j].type = OUTRANGE;
  1026.         else if (!autoscale_ly) {
  1027.             if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
  1028.             cp->points[j].type = OUTRANGE;
  1029.         } else
  1030.             cp->points[j].type = INRANGE;
  1031.         }
  1032.         j++;        /* next valid entry */
  1033.     }
  1034.     /* insert invalid point to separate curves */
  1035.     if (j < cp->p_count) {
  1036.         cp->points[j].type = UNDEFINED;
  1037.         j++;
  1038.     }
  1039.     first_point += num_points;
  1040.     }                /* end while */
  1041.     cp->p_count = j;
  1042.     cp_extend(cp, j);
  1043. }
  1044.